Classical percolation transition 
in the diluted two-dimensional S = 1/2 Heisenberg antiferromagnet 
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The two-dimensional antiferromagnetic 5 = 1/2 Heisenberg model with random site dilution is 
studied using quantum Monte Carlo simulations. Ground state properties of the largest connected 
cluster onlxL lattices, with L up to 64, are calculated at the classical percolation threshold. In 
addition, clusters with a fixed number N c of spins on an infinite lattice at the percolation density are 
studied for N c up to 1024. The disorder averaged sublattice magnetization per spin extrapolates to 
the same non-zero infinite-size value for both types of clusters. Hence, the percolating clusters, which 
are fractal with dimensionality d = 91/48, have antiferromagnetic long-range order. This implies 
that the order-disorder transition driven by site dilution occurs exactly at the percolation threshold 
and that the exponents are classical. The same conclusion is reached for the bond-diluted system. 
The full sublattice magnetization versus site dilution curve is obtained in terms of a decomposition 
into a classical geometrical factor and a factor containing all the effects of quantum fluctuations. The 
spin stiffness is shown to obey the same scaling as the conductivity of a random resistor network. 



PACS numbers: 75.10.Jm, 75.10.Nr, 75.40. Cx, 75.40. Mg 



I. INTRODUCTION 

The two-dimensional (2D) Heisenberg antiferromag- 
net on a square lattice can be driven through a quan- 
tum phase transitior£H3 by, e.g., introducing; frustrating 
interactions^ or by dimerizing the lattice.cl It has also 
been believed that a non-trivial (quantum) phase tran- 
sition could be achieved by diluting the system, i.e., by 
randomly removing either sitesBti or bonds.Erlij The site 
dilution problem is of direct relevance in the context of 
antiferromagnetic, layered cuprates doped with nonmag- 
netic impurities OT13 Diluted Heisenberg models are also 
of more general interest, as systems in which the com- 
bined effects of disorder and quantum fluctuations can be 
studied with a variety of analytical and numerical meth- 
ods. The single impurity problem has been studied ex- 
tensively and is now rather well understood Ji-3 Systems 
with a finite concentration of impurities are much more 
difficult to treat, both analytically and numerically. The 
location and nature of the phase transition driven by di- 
lution is therefore still controversial. 

An early quantum Monte Carlo (QMC) study of the 
temperature dependence of the correlation length gave 
a bound p c > 0.2 for the critical fraction of removed 
sites above which thej-Jong-range order vanishes in the 
2D Heisenberg modelB QMC calculations in the ground 
state indicated p c « 0.35.Q Various analytical treatments 
have given results for p c ranging from 0.07 to 0.30BeI 
These estimates for the critical hole concentratioru-ajrej 
below the classical percolation threshold p* w 0.407,y't£l 
and hence the phase transition would be caused by quan- 
tum fluctuations. A critical hole density much smaller 
than the percolation densitv-was also found in the bond 
diluted Heisenberg model. Q'EEl 

An unusual type of quantum phase transition in t. 
site diluted system was recently claimed by Kato et al. 



They carried out QMC simulations of larger lattices at 
lower temperatures than in previous works and found ev- 
idence of the critical dilution coinciding with the classical 
percolation point; p c — p* . In spite of this, they argued 
that the transition is a non-trivial quantum phase tran- 
sition, which would be a consequence of the fractal clus- 
ters at p* being quantum critical (i.e., with algebraically 
decaying spin-spin correlation function). This leads to 
non-classical critical exponents, which furthermore were 
found to be non-universal, dependent on the spin S of 
the magnetic sites (approaching the classical values when 
S — > oo). Although such behavior violates the standard 
notions of universalis, it cannot be completely excluded 
for random systemslia However, in another recent study 
the spin correlations of the percolating 2D Heisenbepg 
model with 5=1/2 were analyzed in greater detail.ta 
It was confirmed that p c = p*—hixL-. in conflict with 
the quantum criticality scenario JiZIEJ'El strong evidence 
was presented of a transition driven solely by percola- 
tion. The exponents should then be identical to those of 
classical percolation for all S. 

This paper presents details of the QMC studies high- 
lighted in Ref. [l9| and introduces further evidence that 
the order-disorder transition in the diluted 2D Heisen- 
berg model indeed occurs exactly at p* and is classical. 
The stochastic series expansion (SSE) QMC methodE3~Efl 
is used to study the ground state of both site and bond di- 
luted systems at their respective percolation points. Site 
diluted systems are also studied for the whole range of 
hole concentrations p < p* . Particular emphasis is put on 
the importance of carefully controlling potential sources 
of systematic errors in the simulations. In studies of dis- 
ordered systems these issues are much more serious than 
for clean systems, because of the necessity to carry out a 
large number of relatively short simulations for different 
samples (in order to obtain accurate disorder averages). 
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Procedures developed to accelerate the equilibration, and 
to detect possible remaining effects of insufficient equili- 
bration and finite temperature, are discussed here and 
constitute an important part of the paper. 

The main physics questions addressed and results ob- 
tained are summarized as follows. At the percolation 
point, the infinite clusters on .a- 2D lattice have a frac- 
tal dimensionality, d = 91/48.E£l An antiferromagnet at 
this special point could in principle be either classically 
critical (if there is long-range order on the fractal clus- 
ters), quantum critical (with power-law decaying spin- 
spin correlation function on the clusters), or quantum 
disordered (with exponentially decaying correlations on 
the clusters). In the last of these cases, the phase tran- 
sition would occur at a dilution fraction less than the 
percolation density, whereas it coincides with the per- 
colation point in the other two cases. In order to de- 
termine which of the three qualitatively different ground 
states is realized in the percolating cluster of the stan- 
dard Heisenberg model, the sublattice magnetization is 
calculated for the largest cluster on L x L lattices at the 
percolation density, with L up to 64. In addition, clus- 
ters of fixed size N c without boundary imposed shape 
constraints (i.e., on an infinite 2D lattice) are studied 
for N c up to 1024. The sublattice magnetization is av- 
eraged over several thousand samples and extrapolated 
to infinite size. The same non-zero value is obtained for 
both types of clusters, showing consistently that they 
are long-range ordered. Self-averaging is demonstrated 
by studying sample-to-sample distributions of the sub- 
lattice magnetization. The existence of long-range order 
on the percolating clusters implies that the order-disorder 
transition driven by dilution occurs exactly at the perco- 
lation threshold and that the critical exponents are clas- 
sical. The same qualitative behavior is found for site and 
bond dilution, but the long-range order is substantially 
weaker in the bond diluted system. 

In order to reliably calculate the experimentally inter- 
esting sublattice magnetization M as a function of the 
site dilution fraction p for all < p < p* , a decompo- 
sition of M(j>) into a classical and a quantum mechan- 
ical factor is used. The classical factor, which contains 
the singular behavior at p — p* , can be easily evaluated 
by classical Monte Carlo. The critical exponent eavern- 
ing its asymptotic p — > p* form is known exactly.E3 The 
quantum mechanical factor is calculated using QMC sim- 
ulations of the largest cluster on L x L lattices. It is only 
weakly dependent on the dilution. The whole M(p) curve 
is determined to an accuracy of a few percent. 

The spin stiffness is also calculated. Based on known 
results for the classical Heisenberg modeled and the long- 
range order found here in the percolating clusters, it is 
argued that the stiffness should obey the same scaling 
as the conductivity of a random resistor network at and 
close to the percolation threshold. The numerical results 
are fully consistent with the known conductivity expo- 
nent. 

The outline of the rest of the paper is the following: 




FIG. 1. A 64 x 64 lattice randomly diluted at p = p* . The 
solid circles indicate magnetic sites belonging to the largest 
connected cluster (note that periodic boundary conditions are 
applied). The other magnetic sites are shown as open circles. 



In Sec. II the various types of diluted Heisenberg lat- 
tices are defined, and the application of the SSE simula- 
tion algorithm to these systems is discussed. The proce- 
dures developed for controlling potential systematic er- 
rors arising from insufficient equilibration and finite tem- 
perature are also introduced here. Simulation data illus- 
trating the convergence criteria are presented in Sec. III. 
In Sec. IV the sublattice magnetization of percolating 
clusters is studied, both for site and bond diluted sys- 
tems. In Sec. V the full sublattice magnetization versus 
site-dilution curve is calculated. Results for the spin stiff- 
ness are presented in Sec. VI. The paper concludes with 
a summary and discussion in Sec. VII. 



II. MODELS AND METHODS 

The antiferromagnetic 5 = 1/2 Heisenberg model on 
several types of random 2D lattices will be considered. 
In all cases, the Hamiltonian can be written in the form 



H 



b=l 



i(b) ■ »i(6)) 



(1) 



where & is a bond index corresponding to two interact- 
ing nearest-neighbor spins i(b),j(b) and Nf, is the total 
number of bonds on the lattice. On a site diluted lat- 
tice a fraction p of the sites are empty (holes) and the 
rest are occupied by spins. Bonds exist between all occu- 
pied nearest neighbor sites. On a bond diluted lattice all 
sites are occupied and nearest neighbors interact with a 
probability p. Note that a diluted lattice typically con- 
tains isolated (free) spins that are not interacting with 
any other spins. They have to be specified in addition to 
the list of bonds {i(b), j(b)} in the Hamiltonian (fil). 



2 



0.0 1 — = ^ ■ 1 ■ 1 ■ 1 — ^ 1 — 1 

0.0 0.2 0.4 0.6 0.8 1.0 

N/L d 

c 

FIG. 2. Distribution of the size of the largest cluster on pe- 
riodic L x L lattices for L = 16, 32, and 64. The probability 
P(N C ) of cluster size N c is graphed versus N c /L d , showing 
scaling with the fractal dimension d — 91/48. Note the struc- 
ture at N c /L d ~ 0.38, which corresponds to lattices where 
instead of one dominant large cluster there are two of ap- 
proximately half the size. 

A. Diluted lattices 

For lattices with N = L x L sites and periodic bound- 
ary conditions, random magnetic configurations (sam- 
ples) are generated by filling each site with probability 
1 — p. The actual number of magnetic sites is hence 
not fixed, but the fluctuations in the density decrease 
as 1/L. The percolation density p — p* is of special 
relevance. According to the most recent simulation, E2l 
p* = 0.40725379(13). Here the value p* = 0.4072538 will 
be used. The largest cluster of connected magnetic sites 
is of particular interest and its properties will be stud- 
ied separately from those of the full lattice. The number 
of spins belonging to the largest cluster is denoted by 
N c . At p — p* , in the limit L — > oo, this cluster is frac- 
tal, with the fractal dimension d known rigorously to be 
d = 91/48E] For large L the average (N c ) ~ L d , and N c 
is therefore typically considerably smaller than the total 
number of spins on the lattice. One can therefore reach 
larger cluster sizes in the QMC simulations by removing 
the spins that do not belong to the largest cluster. This 
will be done here in order to study the clusters for L as 
large as 64. An example of a a diluted lattice and its 
largest cluster is shown in Fig. [j] 

The largest cluster on a lattice at p — p* exhibits 
strong size fluctuations, as shown in Fig. |^. As an al- 
ternative to approaching the infinite fractal lattice as a 
function of L with fluctuating N c , clusters with fixed N c 
and shapes not restricted by lattice boundaries will also 
be studied. Such clusters are constructed starting from 



FIG. 3. A cluster with N c = 1024 sites constructed on an 
infinite 2D lattice at the percolation density. 

an infinite 2D lattice with only a single filled site. The 
four neighbors of this site are filled at random with prob- 
ability 1— p*. In the next step the neighbors of those sites 
that were filled are in turn filled with probability 1 — p* , 
taking into account that sites that were previously visited 
should not be visited again. This procedure is repeated 
until no new sites can be filled that are connected to the 
cluster, i.e., the nearest neighbors of all sites in the clus- 
ter have already been visited. If the cluster is completed 
before it reaches the desired size N c , or if the size exceeds 
N c , the cluster building is restarted. The process is re- 
peated until a cluster is completed exactly at the size N c . 
This method of constructing fixed-size clusters becomes 
very time consuming for large N c , but it works well for 
sizes N c < 1024 considered here. An example of this type 
of cluster is shown in Fig. |[ 

In the case of .bond dilution, the percolation point is 
exactly p* — 1/2 Jia For L x L lattices this probability can 
be realized for any L and therefore random lattices with 
exactly half of the bonds removed will be considered in 
calculations at the percolation threshold. 

B. Quantum Monte Carlo algorithm 

The r-§SE approach to QMC simulation of lattice 
models£3 has been discussed in detail in previous papers. 
Its application to the Hciscnbcrg model is discussed in, 
e.g., Refs. p3| , p4|j2rj . Its effectiveness for various ordered 
and disordered systep^ has recently been documented 
by several groups .E3c3 Here only a very brief summary 
will be given, in order to facilitate the subsequent discus- 
sion of procedures developed for efficient equilibration 
and ground state convergence for disordered systems. 

In order to apply the SSE method to the Heisenberg 
model, the Hamiltonian ([!]), with J = 1 hereafter, is first 
written as 



H = -J2[Hi, b -H 2 , b ], (2) 

6=1 

where the pair interaction has been divided into terms 

Hl,b = \ — &i(b)Sj(b)i (3) 
H2,b = \ [St( b ) 5*7(6) + S i(b) S t(b) 1 ' ( 4 ) 

which are diagonal and off-diagonal, respectively, in the 
basis {|a)} = {\Sl , S|, . . . , Sff}} used in the simulations. 
A constant has been added to the diagonal part, and as 
a result all non- vanishing matrix elements equal 1/2 and 
correspond to operations on anti-parallel spins. 

The SSE algorithm is based on importance sampling 
of the terms of the partition function Z = Tr{e~@ H } 
written in a truncated Taylor expansion form: 

v^v^ 8 n (M -n)\ , 
Z = E E Ml ("I II H ^M)- (5) 

a Sm i=l 

The summation symbol Sm refers to a sequence of M 
operator- index pairs, 

Sm = [ai,b{\, [02,62], • ■ • ■ [aM,b M ], (6) 

where € {1, 2}, bi E {1, . . . Af,}, corresponding to the 
operators H au \ >i in (||) and ([|), or [ai,bi] — [0,0], cor- 
responding to an identity operator i?o,o = I- This new 
operator has been introduced in order for the summation 
over all Sm i n (r) to imply summation of the Taylor ex- 
pansion of e - ' 3 up to order M. The order of a given 
term corresponds to the number of non-[0, 0] elements in 
Sm, which is denoted by n in (|^). It has been assumed 
that the lattice is bipartite. All the signs arising from 
the off-diagonal operators i?2.& in (§) then cancel in the 
non-vanishing terms of (|E|) and the expansion is hence 
positive-definite. The cut-off M can be easily adjusted 
so that n never reaches M during the simulation. The 
truncation then does not constitute an approximation, 
and SSE simulation results are thus exact to within sta- 
tistical errors. As will be explained further below, M has 
to be chosen proportional to N/3. 

For the sampling of the terms (a, Sm) an efficient algo- 
rithm with three basic updates has been developed. The 
first update involves only diagonal operators. The se- 
quence Sm is scanned from i = 1 to M, and for each 
element [aj,&j] with dj = or at = 1 a substitution 
[0,0] <-> [l,6i] is attempted. The Metropolis acceptance 
probability can be easily calculated from Eq. (g), tak- 
ing into account also that an update in the — > direction 
is allowed only if the spins at the tentative bond bi are 
anti-parallel after operation with the previous i — 1 op- 
erators. An accepted single-operator update changes the 
expansion power n in (|5|) by ±1. 

The second update is a more complicated cluster- 
type update which operates at fixed n and simultane- 
ously changes the operator-type index of several elements 



{i}. The set {i} forms an "operator loop", the size of 
which can be very large. For each i the substitution 
[1, bi] <-> [2,bi] can be carried out without changing the 
configuration weight. The whole sequence Sm can be 
uniquely decomposed into a number of operator-loops, 
which can be updated independently of each other with 
probability 1/2. Details of this operator-loop update are 
discussed in Ref. 

Spins in the state \a) that are not acted upon by any 
operator in Sm are flipped with probability 1/2. Apart 
from isolated spins on a diluted lattice, such free spins 
appear frequently only at high temperatures. 

With the three updates described above — single- 
operator (or diagonal), operator- loop, and spin flip - 
the SSE method is completely grand canonical, i.e., all 
magnetization and winding number sectors are sampled. 
In systems with no isolated spins, the spin flip is strictly 
not needed, but it is still useful at high temperatures. 

The simulation is started with an arbitrary state \a) 
and a short index sequence containing only [0, 0] elements 
(any M will do — in the work discussed here M = A&/4 
was typically used). M is adjusted during the equili- 
bration of the simulation, so that it always exceeds the 
maximum n reached (by, e.g., 20%), and is thereafter 
kept constant. A Monte Carlo step (MC step) consists of 
a full sweep of single-operator updates followed by con- 
struction and updates of all operator loops. After this, 
free spins in the state |a) are flipped with probability 1 /2. 
Further details of the sampling procedures have been de- 
scribed in Refs. |I| and ||(|. 

In the computer, an operator [a, 6] can be represented 
by a single 4-byte integer. In addition, in the cluster 
update four integers are needed to store each operator 
element-iii Sm with their pointers to other elements in 
the listJp The total memory requirement is thus 20 x M 
bytes,Eil plus a few arrays the sizes of which scale lin- 
early with the system size N. The number of operations 
needed for carrying out one MC step scales as M, i.e., is 
proportional to Nj3. 

Observables are typically measured after every MC 
step (it is often practical to do the calculations in com- 
bination with the single-operator update). Estimators 
for various expectation values of interest in the context 
of the Heisenberg model have been discussed in Ref. |2^. 
In the present work, the most important quantity is the 
staggered structure factor, defined on the whole L x L 
lattice as (for a given disorder realization, with S z — 
on the non-magnetic sites) 

S(n,n) = ^^J2(-ip^S! S j ^, (7) 

and on the largest cluster C (or the single cluster on the 
infinite lattice) 



S c (tt,it) = 
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Disorder averaged sublattice magnetizations are defined 
in terms of the structure factors according to 



<m 2 ) = (3S(n,n)/N), 
{ml) = (35 c (7r,7r)/iV c ) 



(9) 
(10) 



where, in the standard way, El the factor 3 accounts for 
rotational invariance in spin space. The order parameter 
m c defined on a cluster will hereafter be referred to as 
the cluster magnetization. 

The spin stiffness will also be discussed. For the non- 
random 2D Heisenberg model with periodic boundary 
conditions it is defined asLJ 



Ps = 77 



3 1 d 2 E {cp) 
2 L 2 dS 2 ' 



(11) 



where is a twist under which the interaction on all 
bonds in one lattice direction is modified according to 



Si ' So 



(12) 



where R is the matrix rotating the 3-component spin vec- 
tor Sj by an angle <f> around the spin-z axis. The stiffness 
can be expressed in terms of the winding number of the 
SSE configurations. The winding number W a , a = x,y, 
is the net number of times spin currents wrap around the 
system in the lattice direction a, i.e., 



W a = (2V* - Nfi/L, 



(13) 



where N£" and are the number of events in the propa- 
gation with the SSE operator sequence Sm in which spin 
is transported to the "left" and "right" along the a di- 
rection. The winding numbers hence take integer values 
0, ±1, ±2 . . .. The stiffness is given by 



(14) 



which can be averaged over the two directions a — x,y. 

For random systems the situation is complicated by 
the fact that the stiffness can vary locally, whereas the 
winding number estimator ( |l4|) is a global quantity char- 
acterizing the rigidity of the system as a whole (i.e., the 
energy increase due to changed boundary conditions). 
This global stiffness is still an important quantity, how- 
ever. One can easily prove that it is equivalent to an 
average stiffness: In a clean system, the definition ( Jl~T| ) 
can clearly be replaced by a definition where the twist 
( |l2| ) is only applied on a single boundary column (which 
has L interacting pairs); 



3 1 d 2 E (<f>) 

2L? d$ 2 ' 



(15) 



The boundary twist here is related to the twist in the 
first definition ( pT| ) by $ = L(f>. If this definition is used 
for a diluted system one still obtains the same expression 
(|l4|) in terms of the squared winding number, regardless 



of which column is taken as the boundary to which $ 
is applied. This is because the spin currents wrapping 
around the system have to go through all L columns. 
The number of interacting pairs on the boundary col- 
umn can depend on which of the L possible columns is 
used, however, and the currents are therefore distributed 
unequally among the bonds although the same net cur- 
rent passes through all columns. This reflects the local 
variation in the rigidity of individual bonds. The stiff- 
ness defined according to the equivalent definitions (|ll]) , 
(|l5|), and (|l4|) is hence the average over all bonds of an 
arbitrary column. In the case that there is no cluster 
wrapping around the system in cither the x or y direc- 
tion, the corresponding winding number is always zero 
and the stiffness in that direction vanishes. Recent dis- 
cussions of the stiffness of disordered quantum systems 
can also be found in Refs. ||J and [35|. 

The bond energy, including the constant added in (j^), 
is obtained in SSE simulations according to 



E b = -(H hb + H 2jb ) = -(n b )//3 ) 



(16) 



where n b is the number of elements [1, b) and [2, b] in Sm- 
Hence, the average expansion power (n) — \E\(3, where 
E is the total internal energy. One can also show that 
the heat capacity C = (n 2 ) — (n) 2 — (n), and hence the 
width of the distribution of n is ~ (n) 1 / 2 at low temper- 
atures. This is the reason why the Taylor expansion can 
be truncated at M ~ N/3. 



C. Convergence Issues 

In QMC studies of random systems, disorder-averaged 
expectation values of the form 



(17) 



normally have to be estimated using only a small subset 
of all Nfj disorder realizations. In addition, the individ- 
ual expectation values (A)i are not evaluated exactly but 
are associated with statistical errors. Typically (A)i is a 
simple operator expectation value [such as the staggered 
structure factor (R) or (||)] which has an estimator that 
is linearly averaged over the importance sampled QMC 
configurations. In principle, the most efficient way to 
estimate the disorder average (0) would then be to gen- 
erate only a single QMC configuration for each randomly 
selected disorder realization, so that each term contains 
both sources of fluctuations (sample-to-sample and QMC 
statistical). The final statistical error can be estimated 
in the standard way using data binning (in order to ap- 
proach a Gaussian distribution from which the standard 
deviation of the average can be calculated). However, in 
practice this approach is not feasible since the simulations 
have to be properly equilibrated for each disorder real- 
ization before the QMC configurations can be used for 
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MC steps 

FIG. 4. Illustration of the /3-doubling scheme used for equi- 
libration and convergence to the ground state. The horizontal 
line segments represent MC steps carried out at the corre- 
sponding inverse temperatures /3 = 2™. No data is collected 
during the N e steps corresponding to the unfilled segments. 
Averages over the two solid segments of length N m are stored 
separately for each /3. 



averages. If a large number of MC steps are needed for 
equilibration it would clearly not be optimal to make use 
of only a single configuration. An accurate estimation of 
the optimum number of configurations would require de- 
tailed knowledge of equilibration times, autocorrelation 
times, and the statistical distributions of the estimators. 
In practice, it is rarely worthwhile to investigate these 
in detail (it would require an effort rivaling that of the 
actual simulations). In any case, the simulations should 
be relatively short so that many disorder realizations can 
be studied. Furthermore, the simulations should not be 
dominated by equilibration. The number of MC steps 
used for sampling expectation values should therefore be 
at least of the same order as the number of steps used 
for equilibration. 

Another important issue is temperature. In order to 
study ground state properties with the SSE method, a 
sufficiently high inverse temperature (3 must be used. 
In diluted systems, especially close to the percolation 
point, different parts of a large cluster may be connected 
only weakly, through essentially one-dimensional narrow 
paths (several examples of which can be seen in Figs. [l| 
and ||). Such "weak links" can lead to correlations that 
develop only at very low temperatures. One can there- 
fore expect that in order to reach the ground state much 
higher (3 values have to be used than for undepleted 2D 
systems. 

Remaining temperature effects and insufficient equi- 
libration are two potential sources of systematic errors 
in the simulations, and these have to be controlled 
very carefully. The following scheme has been devel- 
oped in order to check for both equilibration and tem- 
perature effects: For each disorder realization, simula- 
tions are carried out at inverse temperatures (3 n = 2™, 
n = 0, 1, . . . , n max . Starting with n = {(3 = 1), a 
number N e of MC steps are first carried out for equili- 
bration. Expectation values are sampled during the fol- 
lowing N m = 2N C steps. At the same temperature, N c 
additional steps are carried out during which no expecta- 



tion values are sampled, again followed by N m sampling 
steps. The second segment of N e + N m steps is a direct 
continuation of the first one, so that the effective num- 
ber of equilibration steps for the second sampling seg- 
ment is four times that for the first one. A disagreement 
between the results of the two sampling segments then 
implies that the simulation at the level of the first seg- 
ment is not sufficiently equilibrated, and the second seg- 
ment may also be affected. If the results agree, one can 
conclude that at least the second segment should have 
equilibration errors that are smaller than the statistical 
errors (although this should also be verified by comparing 
simulations with different N m , which will be done below). 
Since the fluctuations of the results of short simulations 
are large, the agreement between the two segments can 
of course be checked only in averages over large numbers 
of simulations of different disorder realizations. 

The /3-doubling scheme is illustrated in Fig. [|. Note 
that simulations at subsequently lower temperatures can 
be started using the last SSE configuration generated at 
the previous temperature. An equilibrated configuration 
at (3 will have an SSE sequence length M approximately 
twice that in the previous run at (3/2. Therefore, in or- 
der to further accelerate the equilibration at low tem- 
peratures, the starting sequence used is the previous Sm 
doubled, i.e., 

S 2 m = [ai,h],..., [a M ,b M ][a M ,b M ], . .., [ai,b x \. (18) 

Especially at low temperatures, where the system is al- 
most converged to the ground state, the doubled SSE 
configuration should be very nearly distributed according 
to the equilibrium distribution at the new (3. With the 
reversed order of the second set of M operators in jl^), 
the initial 5*2 m pAfvays has zero winding number, which 
can be expected^ to be a slightly better starting point 
than the alternative one with twice the winding number 
(in practice, the difference in performance is minor). 

Expectation values calculated for all n max + 1 values 
of (3 are stored on disk, so that the convergence to the 
ground state can be checked. Ideally, the number of (3- 
doublings should be large enough that there are no sta- 
tistically significant differences between the results for 
fimax = 2™ max and (3 = 2™ max_1 . Since the asymptotic 
convergence is exponential, the results at (3 ma x should 
then have no detectable temperature effects at the level 
of the statistical errors. 



III. CONVERGENCE TESTS 

In this section test results for equilibration and ground 
state convergence according to the /3-doubling procedures 
described in the previous section are presented. Dilu- 
tion fractions close to the percolation point can be ex- 
pected to be the worst with respect to slow (3 conver- 
gence. This is because for p < p* the largest clusters are 
two-dimensional and more compact than at p* (i.e., they 
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FIG. 5. Test results for the convergence of the sublattice 
magnetization of the largest cluster on 32 x 32 lattices at p* . 
The number of MC steps for averaging data for each point 
was N m = 2. The open and solid circles correspond to the 
first and second data collection segment, respectively. The 
results are averages over 10 4 disorder realizations. 



have less "weak links"), and for p > p* the cluster size 
does not diverge with L. Site diluted systems exactly at 
the percolation point are considered here. 
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FIG. 6. Dependence of the calculated cluster magnetiza- 
tion on the number of MC steps in the data collection seg- 
ments, N m = 2 fe , at two different inverse temperatures (re- 
sults for L = 32 lattices at p* , averaged over 10 4 samples). 
The open and solid circles correspond to the first and second 
data collection segment, respectively. 



A. Equilibration 

The equilibration of the simulations will first be illus- 
trated by results for L = 32 systems obtained with differ- 
ent N e and N m = 2N e . Fig. || shows results for the dis- 
order averaged cluster magnetization when the segments 
are very short; N e — 1 and N m = 2. At the highest 
temperature, (3=1, the two segments give results that 
agree within statistical errors, but as the temperature is 
lowered the results begin to differ considerably. At still 
lower temperatures the results again converge and be- 
come statistically indistinguishable at j3 = 1024 in this 
case. The good agreement here can be explained by the 
fact that low-temperature simulations in the /3-doubling 
procedure start from configurations which already have 
a rather long history at higher temperatures, which in 
combination with the trick of doubling the SSE operator 
sequence produces almost equilibrated initial configura- 
tions when the system is nearly in its ground state. 

Fig. |6| shows how results at an intermediate and low 
temperature depend on the number of MC steps in the 
data collection segments. At (3 = 32, the first data seg- 
ment converges after N m > 16, whereas the second seg- 
ment appears to be converged already for N m = 4. At 
/3 = 2048, the results for the two segments agree statis- 
tically for all N m , and the averages show no discernible 
dependence on N m . Hence, an agreement between the 



two segments indeed appears to be a good indication 
of sufficient equilibration. Since the convergence is the 
slowest at intermediate temperatures, a very safe conser- 
vative check of low-temperature equilibration should be 
that the two segments agree at all temperatures. For the 
final result, the segments can then be averaged in order to 
improve the statistics. However, this typically leads only 
to a modest reduction of the error bars (i.e., significantly 
less than the reduction by y/2 expected for independent 
data) since the statistical errors are dominated by fluc- 
tuations between the disorder realizations. The fact that 
sample-to-sample fluctuations dominate can also be seen 
clearly in Fig. ||, where the error bars decrease much 
slower than by \[2 for successively higher k. 

The results presented here indicate that even ex- 
tremely short simulations give results void of non- 
equilibration effects at low temperatures. However, 
longer runs were used to produce some of the data pre- 
sented in this paper. The main reason for this is that al- 
though unbiased disorder averages of the form (|l7|) can be 
obtained with short simulations, large statistical errors 
in the individual expectation values can be problematic 
when considering non-linear functions of the expectation 
values (such as their typical values) or their complete sta- 
tistical distributions. One then has to demand that the 
statistical errors of the individual expectation values are 
much smaller than the width of the distribution of the 
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FIG. 7. Distribution of the cluster magnetization of 32 x 32 
lattices at p* for different lengths of the data collection seg- 
ments. The inverse temperature /3 = 4096, and 10 disorder 
realizations were used for each N m . 

exact expectation values. An example of how statistical 
errors can distort distributions is shown in Fig. ^ where 
histograms of the cluster magnetization are compared for 
six different simulation lengths. Both the data collec- 
tion segments were used for calculating the individual 
expectation values, i.e., the number of measurements for 
each realization is 2N m . The histograms become signifi- 
cantly narrower as the number of MC steps is increased. 
The distribution is not completely converged even for 
the longest simulation considered here (N m — 64), but 
the relatively small differences between N m = 32 and 64 
suggest that the N m = 64 result is close to the exact dis- 
tribution. Note that the first moment of the distribution, 
i.e., the linear disorder average ([l7|), is the same within 
statistical errors for all N rn (which was demonstrated at 
(3 = 2048 in Fig. g). 

In the calculations discussed in the following sections, 
N m between 100 and 250 was used to ensure that reli- 
able distributions could be obtained at the percolation 
point. For p < p* , where the full distributions are not as 
important, N m = 50 was typically used. Since effects of 
insufficient equilibration are undetectable even in much 
shorter simulations the results should definitely be void 
of any bias of this nature. 

B. Ground state convergence 

Already the results shown in Fig. g demonstrate that 
very low temperatures are required in order to converge 
the sublattice magnetization to its ground state value. 
For L = 32, a /3-value higher than 2000 is needed to 
eliminate temperature effects within the statistical er- 
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FIG. 8. Cluster magnetization ratios vs inverse tempera- 
ture for L x L lattices. The number of samples used for 
averaging the results was 88000, 21000, 9000, and 2500, for 
L = 8, 16, 32, and 64, respectively. 

rors. In order to more accurately study remaining ef- 
fects at low temperatures it is useful to calculate ratios 
(m^(Pi)) I (ml((3j)) of the squared cluster magnetization 
at different temperatures. The relative statistical errors 
are smaller in the ratios than in the absolute values, since 
the sample-to-sample fluctuations cancel when the same 
realizations are used at all temperatures. Fig. || shows 
results for systems with L = 8, 16, 32 and 64, which were 
simulated with (3 up to /3 max = 256 x L. The ratios, 
with the data at the respective /3 max in the-denominator, 
were analyzed using the bootstrap methodED in order to 
obtain accurate estimates of the error bars. For L = 8 
and L = 16, the results at /3 max and /3 max /2 do not dif- 
fer within statistical errors and hence the result at /3 max 
should not have any temperature effects left at this preci- 
sion level. The L — 32 and 64 results are not completely 
converged to the ground state, however. The exponen- 
tial low-temperature convergence seen for all the system 
sizes indicates that the remaining temperature effects at 
/3 max should only lead to an error that is smaller than 
the difference between the ratios at /3 max and /3 max /2. 
Hence, the under-estimation of the sublattice magnetiza- 
tion should be less than 0.2% for L = 32 and less than 
0.5% for L = 64. These upper bounds for the systematic 
errors are of the same magnitude as the respective statis- 
tical errors in {m 2 c ) (which unlike the ratios also include 
contributions from sample-to-sample fluctuations). The 
remaining small temperature effects should therefore not 
substantially affect the finite-size scaling of the sublattice 
magnetization (to be discussed in the next section). 

Fig. U shows magnetization ratios for fixed-7V c clus- 
ters, with max = 32 x N c . In this case there are small 
but detectable differences between the results at /3 max 
and /3 max /2 for all system sizes, except for N c = 1024 
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where the statistical error is larger than the difference. 
Again, the maximum relative systematic errors remain- 
ing at pmax are similar in magnitude to the statistical 
fluctuations in (m 2 ) and can only have very minor ef- 
fects on the finite-size scaling. 

The j3 needed for ground state convergence decreases 
rapidly away from the percolation point, and therefore 
the p < p* results for L x L lattices discussed in Sees. V 
and VI are completely converged even for L = 64. 



IV. LONG-RANGE ORDER IN PERCOLATING 
CLUSTERS 

In this section the ground state sublattice magnetiza- 
tion of the percolating cluster is investigated in detail. If 
it remains finite in the thermodynamic limit, the order- 
disorder transition driven by dilution must necessarily 
occur exactly at the classical percolation density. To see 
this, consider the sublattice magnetization (||) of the di- 
luted L x L lattice. Its disorder average can be written 
as a sum of contributions from all the clusters k on the 
lattices as, 
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FIG. 9. Magnetization ratios vs inverse temperature for 
fixed-iV c clusters. The number of samples was 43000, 15000, 
10000, 3000, and 1100, for iV c = 64,128,256,512, and 1024, 
respectively. 



\ k I 



(19) 



In the thermodynamic limit, only infinite clusters con- 
tribute to this sum, and therefore one only needs to con- 
sider the behavior of the cluster magnetizations m 2 for 
large clusters. If there is long-range order, it is natu- 
ral to assume that the sublattice magnetization is self- 
averaging (a fact that will also be demonstrated explicitly 
below). The individual to 2 values can then be replaced 
by the infinite-size extrapolated average for the largest 



cluster, i.e. 



(m 2 ) 



}, which gives 



N 2 



E^ 2 



(L -> oo). 



(20) 



This expression is identical to the order parameter of a 
classical diluted system, up to the factor (m 2 ) which is 
reduced by quantum fluctuations from its classical value 
1/4 (for an Ising model with S? — ±1/2). If (m 2 ) remains 
finite at p = p* [which is the condition for (EQ) to remain 
valid for all p < p*] the only singular behavior is in the 
classical expectation value and hence the critical behavior 
is that of classical percolation. 

In general, Eq. (^0|) holds for any dilution fraction 
p < p c , where p c in principle may be less than p*. The 
cluster magnetization (m 2 ) 1 / 2 will here be determined at 
percolation, where the infinite clusters are fractal. Dilu- 
tions less than the percolation density, where the infinite 
clusters are two-dimensional, will be studied in the next 
section. 



A. Sublattice magnetization in site-diluted systems 

At the percolation point, the average number of spins 
in the largest cluster on a diluted L x L lattice scales 
asymptotically as (N c ) ~ L d , with d the fractal dimen- 
sion 91/48.E.3 As can be seen in Fig. || the full distribu- 
tion of the size of the largest cluster also scales as L d , 
i.e., the distribution width also diverges as L — > oo. This 
is in sharp contrast to the situation below the percola- 
tion threshold where the size distribution approaches a 
(^-function at a size ~ L 2 . Note, however, that the scaled 
distribution at p* has sharp cut-offs both at the lower and 
upper edge, meaning that also the smallest and largest 
clusters grow as L d . Hence, finite-size scaling of (m 2 ) 
calculated on such fluctuating- N c clusters as a function 
of L is a well defined procedure for extracting the sub- 
lattice magnetization of the infinite fractal cluster. Nev- 
ertheless, the alternative way of approaching the ther- 
modynamic limit with fixed-iVc clusters on the infinite 
lattice is also considered here. An agreement between 
the two calculations will provide additional support to 
the argumentt2l that the percolating cluster is ordered. 

In the rpure 2D Heisenberg model the leading size- 
correctiorH to m 2 Js^ iV -1 / 2 , which can be seen clearly 
in numerical dataEj'Ej In analogy with this, as a scaling 
hypothesis at percolation, the following leading size cor- 
rections are tested here for the largest cluster on L x L 
lattices and fixed- N c clusters, respectively: 



(m 2 ) x +aL- d '\ 
(mDx^imD^ + bN- 1 / 2 . 



(m 2 ) L 



(21) 
(22) 



Fig. [lO] shows results for L up to 64 and N c up to 1024. 
The data are fully consistent with the scaling ansatz, al- 
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FIG. 10. Finite-size scaling of the disorder averaged 
squared cluster magnetization. The results for the largest 
cluster on L x L lattices (solid circles) are graphed vs L~ d ^ 2 
and those for fixed- N c clusters (open circles) vs iVc -1 '' 2 . Sta- 
tistical errors are much smaller than the symbols. The curves 
are cubic polynomial fits. 



though in order to fit all the points a polynomial cubic 
in L~ d l 2 has to be used in both cases (a cubic polyno- 
mial is needed also to .fit high-accuracy data for the clean 
2D Heisenberg mode£3). The infinite-size extrapolated 
values for {m 2 c ) from the two fits agree very well (within 
statistical errors) . The sublattice magnetization is in fact 
quite large, (m c ) = 0.150(2), which is almost precisftbx 
half of the value m = 0.307 for the clean 2D system£3'Ea 

It should be stressed that it is not critical whether or 
not the scaling ansatz assumed here to carry out the ex- 
trapolation of the sublattice magnetization is strictly cor- 
rect or not. Unless the behavior would change dramati- 
cally for even larger systems, a slightly different finite-size 
correction would not significantly affect the extrapolated 
(m 2 ). One could of course argue that a cross-over to a 
qualitatively different -behavior cannot be excluded, as 
indeed has been done.EHl No plausible physical reason for 
such a cross-over has been presented, however. With two 
different boundary conditions for the clusters giving the 
same result for the infinite-size extrapolated sublattice 
magnetization, the most natural scenario must be that 
the percolating cluster is ordered. 

In a disordered system the order parameter is not con- 
stant over the whole system, but depends locally on the 
structure of the lattice. One would, however, expect 
self-averaging, i.e., the sublattice magnetization averaged 
over different regions of an infinite cluster should be the 
same when the size of the regions is sufficiently large. In 
finite systems, self-averaging can be seen in the statisti- 
cal distributions of the individual cluster magnetizations. 
Fig. [ll] shows results at the percolation point for several 
L x L and fixed-iV c systems. As discussed in Sec. Ill, 
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FIG. 11. Distribution of the cluster magnetization of L x L 
lattices (top panel) and fixed- N c clusters on the infinite lattice 
(bottom panel). 



0.100 
0.090 
0.080 
0.070 
0.060 
0.050 
0.040 
0.030 
0.050 
0.045 
0.040 
0.035 
0.030 
0.025 
0.020 
0.015 



N=128 




w-, .• •• • . • 

§^i*v ".-'V ••• 



10 



N=1024 



15 



20 



25 



30 



FIG. 12. Individual squared magnetizations vs the radius 
of gyration for clusters on the infinite lattice. 
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the histograms can be expected to be slightly broadened 
by the statistical fluctuations in the SSE results for the 
individual m;: values. Such effects should, however, be 
minor when the simulations are as long as those used 
for the data shown here (N m = 100 for the L x L lat- 
tices and 250 for the fixed- N c clusters). The widths of 
both types of distributions decrease with increasing sys- 
tem size, which is consistent with vanishing fluctuations 
in the thermodynamic limit. It can also be noted that 
the distributions become more symmetric for larger sys- 
tems — the weak tails visible at the high-end of the dis- 
tributions for small clusters vanish as the system grows. 
The behavior is hence fully consistent with the ^-function 
distribution expected for a self-averaging quantity in the 
thermodynamic limit. 

It is also interesting to study how the cluster magne- 
tization depends on the shape of the cluster. A compact 
cluster is likely to have a stronger order than one which 
has many narrow paths. A natural length scale charac- 
terizing the over-all density of the unconstrained fixed- N c 
clusters is the radius of gyration, 

(\ 1/2 
c i=l;=l / 

where Xi,yi are the (integer) coordinates of the magnetic 
sites. Fig. [l2] shows scatter plots of the cluster magneti- 
zation versus R for two cluster sizes. For N c = 128, one 
can see that the most compact clusters, i.e., those with 
the smallest R, indeed have the largest magnetizations. 
After an initial rapid decrease with R for the smallest 
R, the average magnetization only decreases slowly with 
increasing R, however. The N c = 1024 clusters show a 
similar behavior. There are of course in principle clus- 
ters with very large R that should have much smaller 
magnetizations, but these clusters lack statistical signif- 
icance. The weak independence for the statistically sig- 
nificant clusters is another manifestation of a strongly 
self-averaging sublattice magnetization. 

B. Sublattice magnetization in bond-diluted systems 

For the bond-diluted system only simulations of L x L 
lattices were carried out. Fig. |l3| shows the results for 
the cluster magnetization at the bond percolation point, 
P* = 1/2, graphed in the same way as for the site diluted 
systems in Fig. Eo. Also in this case the scaling to a finite 
sublattice magnetization is evident, but the value of the 
order parameter is smaller than in the site diluted system; 
(m c ) — 0.088(2). The difference can be explained by 
the different local structures of the two types of lattices. 
Although the fractal dimension d-oi the cluster is the 
same for site and bond percolation,Ej the average number 
of bonds per spin is smaller in the bond diluted case 
- 1.121 versus 1.259. This leads to stronger quantum 
fluctuations in the bond-diluted system. The infinite-size 
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FIG. 13. Finite-size scaling of the disorder averaged 
squared cluster magnetization of the bond diluted system at 
the percolation density. The curve is a cubic polynomial fit. 

energy per bond (which reflects the tendency to nearest- 
neighbor singlet formation) is -0.3890(1) and -0.4068(2) 
for site and bond dilution, respectively. 

It can also be noted that for a given L the average 
largest cluster on the bond diluted lattice is w 45% larger 
than on the site diluted lattice. The stronger quantum 
fluctuations and the larger cluster sizes imply that for 
given L a lower temperature has to be used to converge 
the bond diluted system to the ground state. For the 
largest size studied in this case, L = 32, an inverse tem- 
perature (5 = 32768 was used. 

C. Scaling of the full staggered structure factor 

The previous,claims of quantum criticality at the per- 
colation point! 1 H 21 ! were primarily based on a finite-size 
scaling analysis of the staggered structure factor. A log- 
log plot of S(ir,ir) calculated using SSE simulations in- 
cluding all the spins of diluted L x L lattices is shown in 
Fig. The numerical values agree well with those of 
Ref. jr7| . One can, however, expect a barely discernible 
finite-T reduction in the previous L — 48 results because 
the temperature used {(3 = 1000) was not sufficiently low 
for complete converge to the ground state (see Fig. ^ and 
a related discussion in Ref. |l9|) . 

The scaling seen in Fig. [14| is different from the ex- 
pected classical percolation behavior. Given the results 
presented above for the scaling of the cluster magnetiza- 
tion, the deviation from classical behavior for this range 
of system sizes is not surprising, however. Classically, 
the finite-size scaling of S(n, n) is solely the result of the 
divergence of the size of the connected clusters with L. 
In the quantum mechanical case, there is also a factor, 
the sublattice magnetization of the cluster, multiplying 
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FIG. 14. Finite-size scaling of the disorder averaged stag- 
gered structure factor of the full site diluted L x L lattice. 
The line has slope 43/24, expected for classical percolation. 



each cluster size, i.e. 



Nlml 



(24) 



The size-dependence of the average m\ was shown in 
Fig. [l(]. From these results it is clear that there is an 
effect that partially compensates for the growth of the 
cluster sizes Nk in Eq. (|24|), namely, m 2 . decreases with 
increasing cluster size. Hence, for systems where the rela- 
tive size corrections to the cluster magnetization are still 
significant, as is the case for all sizes that can currently be 
reached in numerical simulations, the growth of 5(77, it) 
with L will be slower than for a classical system. This 
explains the slow convergence towards the classical be- 
havior that can be seen in Fig. 
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It can be noted that 
the largest cluster completely dominates the staggered 
structure factor and the curve shown in Fig. |lj changes 
only very little if only the largest cluster is included, i.e., 
S c (tt,it)^ S (it, n)H 



V. DILUTION DEPENDENCE OF THE 
SUBLATTICE MAGNETIZATION 

The doping dependence of the sublattice magne- 
tization of antiferromagnetic cuprates can be mea- 
sured experimentally using NQR, /i + SR and neutron 
scattering. tilll^l Results for the Heisenberg model were 
recently obtained using an improved spin-wave the- 
ory which, however, breaks down close to the perco- 
lation threshold (the critical point is unphysical, lo- 
cated at a bo Jc concentration higher than the percolation 
density) .Bcil Previous QMC calculations of the doping 
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FIG. 15. Size dependence of the classical magnetization for 
systems diluted at one half percent less than the percolation 
density (p = p* — 0.005). The open circles correspond to 
the full cluster sum, Eq. (|^). The solid circles are from the 
average including only the largest cluster, Eq. (pi|). 

dependence were based on Eq. (§).0 Use of this formula 
becomes very difficult close to the percolation threshold, 
however, since the smallness of the sublattice magnetiza- 
tion there is associated with a slow convergence to the 
asymptotic regime in which finite-size scaling is reliable. 
Here a different approach will be taken, based on the 
fact that the sublattice magnetization can be decomposed 
into classical and quantum mechanical factors, which can 
be evaluated separately. This decomposition was already 
discussed in Sec. IV and was written down as Eq. (|2"c|). 
Here the notation M = (to 2 ) 1 / 2 will be used for the dis- 
order averaged sublattice magnetization. Eq. ( pp[ ) can 
then be written as 



M(p) = M qm (p)M c i(p), (L -> oo), 
where M qm is the quantum mechanical factor 



M„ 



and M c i is the classical (geometrical) factor 

\ 1/2 



E^ 2 



(25) 



(26) 



(27) 



In the ordered regime, < p < p* , only the largest cluster 
contributes to this sum in the thermodynamic limit. The 
classical factor can therefore also be obtained as 



M cl = (N c )/N. 



(28) 



Fig. [D] shows the size convergence using the two defi- 
nitions of the classical factor when the dilution fraction 
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FIG. 16. Upper panel: Classical magnetization vs dilution. 
The dashed line shows the small-p form M c ; = 1 — p. Lower 
panel: The magnetization divided by (p* — p) 5//36 graphed vs 
p* — p. The curve is a polynomial fit, with parameters given 
inEq. ©. 



p = p* — 0.005. The single-cluster average (f28|) clearly 
converges faster. It is apparent that a reliable extraction 
of the quantum mechanical sublattice magnetization M 
using the structure factor formula (^|) would be impossi- 
ble in this case, since not even the classical magnetiza- 
tion is in the asymptotic scaling regime for the range of 
system sizes where QMC simulations can be carried out 
(L < 100). The quantum mechanical M can be expected 
to have an even worse scaling behavior, due to effects 
similar to those found for the staggered structure fac- 
tor in Sec. IV-C. The quantum mechanical factor M qm 
can be calculated based on much smaller system sizes, 
however. It was evaluated in the extreme case p = p* 
in Sec. IV, and even there it is as large as 50% of the 
value in the other extreme, i.e., the non-diluted system 
(p = 0). Hence, M qm is only weakly dependent on the 
dilution fraction, and most of the p dependence of M, 
including the singular behavior at p*, is in the classical 
factor M c i. 

The classical magnetization is known to vanish at the 
percolation threshold with the exponent 5 / 36 Jl3 i.e., 



M cl (p^p*)=A cl (p*~p) 5 ^. 



(29) 



suit Md = 1 — p. Numerical values for < p < p* — 0.002 
were obtained here by simulations of lattice sizes as large 
as L — 4096, using the single-cluster estimator (|2q). 
In the example illustrated in Fig. [l^, the results for 
the three largest sizes are 0.43640(4), 0.43636(3) , and 
0.43634(2) (for L = 512, 1024, and 2048, respectively) 
and the L = 2048 result (which is based on 3 x 10 5 
samples) can hence be taken as the infinite-size value 
of M c i(p* — 0.005). Closer to the percolation threshold 
L = 4096 lattices were used. Fig. [16] shows results for the 
whole dilution range. The linear small-p form describes 
the data well for p up to « 0.2. The asymptotic form 
( p9| ) is well reproduced for p* — p < 0.02, with the con- 
stant A c i 0.91. In order to have an analytic expression 
describing M c \ in a wider region around p* , higher-order 
terms can be added to A c i. The following forms will be 
used in combination with fits to the quantum mechanical 
factor in order to obtain expressions for M both close to 
p = and p = p*: 



M cl (p < 0.2) = 1 - p, 
M c i{p* -P < 0.2) = [0.9102 + 3.053(p* -p) 2 

-5.642(p*-p) 3 )](p* -p) 5 / 36 . 



(30) 
(31) 



Note that it is not claimed here that the higher-order 
terms in the form (|3l]) are the correct subleading terms 
of the critical percolation behavior — the purpose is just 
to have an expression that describes the data well in prac- 
tice, within the stated region. 

The quantum mechanical factor can be calculated in 
the same way as was already explained in the case of 
p = p* in Sec. IV, i.e., using the SSE method for the 
largest cluster of connected magnetic sites on L x L lat- 
tices. Here L up to 64 was used for dilution fractions 
p = p* - 6, with S = 0.05, 0.10, . . . , 0.35, and 0.38. SSE 
simulations at p = were previously carried out for L 
up to 16,E3 and were here extended up to L = 64. Since 
the ground state convergence occurs at lower (3 as 5 is in- 
creased, the simulations are faster than at p* and a larger 
number of samples could therefore be studied. The nu- 
ber of samples was > 10 4 for L = 64 and up to 10 6 for 
smaller lattices. The results were extrapolated to infinite 
size using a leading .correction ~ 1/L (as in the case of 
the clean 2D systemE3). The resulting M qm (p) is shown 
in Fig. |l7|, along with two quadratic fits which describe 
the data very well over quite wide ranges of p. The fitted 
forms are 

M qm {p < 0.25) = 0.3072 - 0.134 ■ p - 0.51 • p 2 , (32) 
M qm (p* -p< 0.25) = 0.151 + 0.721 • (p* - p) 

-0.93- (p* -p) 2 . (33) 

The final result for the sublattice magnetization of the 
site diluted Heisenberg model is shown in Fig. [l^. The 
solid circles were obtained by interpolating the numerical 
results for M c i and M qm and multiplying them according 



In the weak dilution limit, one can easily obtain the re- to Eq. (g5j). Forms describing the results well in quite 
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FIG. 17. The quantum mechanical factor (cluster magne- 
tization) vs dilution. The curves are quadratic fits (solid for 
small p and dashed for small p* — p) with parameters given 
in Eq. (B2) and (B3h, respectively. 
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FIG. 18. Sublattice magnetization vs site-dilution (solid 
circles). The curves are parametrized forms discussed in the 
text. The inset shows the behavior close to the percolation 
threshold on a more detailed scale. 



wide regions p < 0.25 and p* — p < 0.25 can be obtained 
by multiplying the corresponding expressions (pOf), ( [32] ) 
and (|3l|), (|33|). The resulting curves are also shown in 
the figure. The initial linear reduction M(p)/M(0) = 
1 — Bop, where the coefficient and its estimated error 
is Bo = 1.44 ± 0.05. At the percolation threshold the 
leading behavior is M(p) — A*(p* — p) 5 / 36 with A* = 
0.137 ± 0.002. 

In Ref. the sublattice magnetization normalized by 
the number of magnetic sites, M'(p) — M(p)/(1 — p), 
which is equivalent to the quantum mechanical factor 
-^qm for small p, was calculated using spin wave the- 
ory with a T-matrix approximation. The initial lin- 
ear weak-dilution form M'(p)/M'(Q) = 1 — B' p, with 
B' Q = 0.691 ± 0.005, was found in that approximation. 
The results obtained here for M(p) correspond to a 
slightly smaller coefficient; B' = 0.44 ± 0.05. Despite 
this difference in initial slope, the relative agreement be- 
tween the full sublattice magnetization shown Fig. [l8| 
and the corresponding spin-wave result is very good up 
to p pa 0.15 (not shown here — see Fig. 10 of Ref. 41). 
For higher p, the spin-wave result falls significantly be- 
low the QMC result until very close to the percolation 
point, where the actual M{p) approaches zero but the 
spin-wave result remains finite.cil For p < 0.35, the re- 
sults shown in Fig. [l8| agxee well with those presented 
previously by Kato et al.Ui Their extrapolations closer 
to the percolation threshold are not reliable, however, 
since they fitted a different, non-classical exponent to de- 
scribe the p — > p* behavior. The estimated accuracy for 
the M(p) curve obtained here is better than 2% over the 
whole range of dilutions (significantly better forp < 0.1). 



VI. SPIN STIFFNESS 

As discussed in Sec. II-B, the spin stiffness p s can 
be obtained in SSE simulations in terms of the wind- 
ing number fluctuations, Eq. (|l4|). At dilution fractions 
p > p* there are no clusters wrapping around the peri- 
odic L x L lattice for large L, and therefore the stiffness 
vanishes identically. Exactly at the percolation thresh- 
old the wrappingrprobability in a given direction is ap- 
proximately 0.52,li3 and the stiffness can then be non- 
zero for finite L. It should vanish in the thermodynamic 
limit, however. For p < p* the wrapping probability ap- 
proaches 1 as L — > oo and in view of the existence of 
antifcrromagnetic long-range order the stiffness can then 
be expected to be non-zero. 

For the classical diluted Heisenberg model the botav- 
ior of the stiffness (or helicity modulus) is known.cj It 
scales in the same way as the conductivity of a ran- 
dom resistot-network, with the conductivity exponent t of 
percolationJl3 According to recent-simulations, the value 
of this exponent is t = 1.310(1) J13 With the long-range 
order present in the percolating clusters of the quantum 
Heisenberg model, as demonstrated in Sec. IV, one can 
expect that the stiffness should behave as in the classical 
model (in analogy with the "renormalized classical" be- 
havior of the clean 2D Heisenberg modelj^l and in view of 
general symmetry arguments). Scaling with the conduc- 
tivity exponent will therefore be tested here. It can be 
noted that the elastic moduli of a diluted classical elastic 
lattice also obey scaling with the conductivity exponent, 
if the force constants are isotropics With non-isotropic 
forces other scalings have been shown to be possible, and 
the critical dilution fraction above which the rigidity, van- 
ishes can in fact be below the percolation density.c3 
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FIG. 19. Spin stiffness of site diluted L x L systems multi- 
plied by I}' v . The curve is a quadratic polynomial in L~ d ^ 2 . 



The elastic moduli of classical percolating lattices 
been studied extensively using numerical methods .c3'l£I 
The primary scaling technique used there will be em- 
ployed here as well. If the analogy with the random re- 
sistor network holds, the disorder averaged_stiffness at 
the percolation point should scale as L -t / 1/ £3 where v is 
the correlation length exponent of percolation, which is 
known exactly; v = 4/3.113 Fig. 19 shows the stiffness of 
the site-diluted Heisenberg model at the percolation dep.. 
sity, multiplied by L l l v , where t/v = 0.9826 was used.c^l 
The data extrapolate to a finite value as L — > oo, and 
hence the results are indeed consistent with the conduc- 
tivity scaling. The average stiffness shown here was cal- 
culated by including only non-zero values of the stiffness 
in a given lattice direction, and was averaged also over the 
two equivalent directions. The simulations included only 
the largest cluster in the system. The fraction of non-zero 
stiffnesses is approximately 0.52 for all L, in agreement 
with the knownMI wrapping probability of clusters in pe- 
riodic systems. 

Fig. |2^ shows the full dilution dependence of the stiff- 
ness extrapolated to infinite system size. The behavior 
is almost linear up to p w 0.15. The fitted linear form 
Ps(p) = 0.1808 — 0.62 • p is shown in the figure, along 
with an analytical result containing terms up to p 2 ob- 
tained using a non-linear er-model and percolation theory 
in Ref. |8|. The a-model approach gives a quantum critical 
point below the percolation point, and the initial fall-off 
is also faster than what is seen in the QMC data. The 
QMC data close to the percolation point are not well 
described by the random resistor network exponent, al- 
though the results at the percolation point indicate that 
this should be the asymptotic form as p — > p* . The criti- 
cal region may be very small, however, making it difficult 
to observe in direct calculations for p < p*. 
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FIG. 20. Infinite-size extrapolated spin stiffness vs dilution 
fraction (circles). The solid line is p a (p) = 0.1808 — 0.62 - p. 
The dotted curve is the scaling form p s {p* — p) ~ (g* — p) 1 , 
and the dashed curve is a non- linear cr-model result .□ 



VII. SUMMARY AND DISCUSSION 

This paper has presented a variety of quantum Monte 
Carlo results showing that the order-disorder transition 
in the diluted S =1/2 Heisenberg model is solely driven 
by classical percolation. This is a consequence of the frac- 
tal clusters at the percolation point having long-range 
antiferromagnetic order. The presence of this long-range 
order was demonstrated by studying the largest cluster 
on L x L lattices, as well as clusters of fixed size N c on the 
infinite 2D lattice. For the infinite-size extrapolated sub- 
lattice magnetization, the same non-zero value was ob- 
tained for both types of clusters. An accurate calculation 
of the sublattice magnetization M versus site dilution p 
was made possible by taking advantage of a factoring 
into classical and quantum mechanical functions, which 
were evaluated separately. The classical factor is identi- 
cal to the magnetization of a classical ferromagnet, and 
was calculated to high accuracy using lattice sizes L up 
to 4096. It contains the critical behavior of M(p) close to 
the percolation point p* . The quantum mechanical fac- 
tor is equivalent to the sublattice magnetization of the 
largest cluster on L x L lattices in the limit L — > oo. It 
remains non-singular as p — > p* and can be reliably ex- 
trapolated using relatively small system sizes (here using 
L < 64). Its infinite-size value grows only by a factor 
2 between p = p* and p — 0. Approximate analytical 
forms describing the numerical M(p) for all p were also 
constructed. The spin stiffness at the percolation point 
was shown to obey the same scaling as the conductivity 
of a random resistor network. 

The conclusions reached in this paper differ from the 
non-universal quantum criticality scpnario which has 
been elucidated in several recent papers £3'E3Eil Accord- 
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ing to this scenario, the fractal clusters at the perco- 
lation point have power-law decaying spin-spin corre- 
lations, which implies that the scaling exponents differ 
from those of classical percolation. It was argued that the 
exponents depend on the spin S of the magnetic sites, sa 
that classical percolation is recovered only for S — ► oo.tll 
Several types of scaling analyses, ha v.e been presented in 
support of this unusual behavior It can be noted, 
however, that only a very small number of systepa. sizes 
were used. Temperature effects, although smalljij may 
also have contributed to making the scaling appear better 
than it would be for real T = data. 

The most serious problem with the finite-size scalings 
carried out in Refs. [17|p0| , and |l] is that even if the per- 
colating clusters are ordered, as they in fact are, the 
staggered structure factor cannot be expected to show 
the asymptotic classical scaling behavior for the range of 
system sizes used (as shown in Fig. [[J]). This is due to 
the strong size dependence of the sublattice magnetiza- 
tion of the clusters (in contrast to the classical case, in 
which the cluster magnetization takes its maximum value 
at T = for any system size). The classical scaling form 
is valid only for systems sufficiently large for the relative 
size corrections to be small, which is the case only for sys- 
tems much larger than those that are currently accessible 
to quantum Monte Carlo simulation. This problem was 
circumvented here by focusing on the cluster-size nor- 
malized sublattice magnetization of the largest cluster of 
the lattice, which in combination with the known scaling 
of the cluster sizes completely determines the asymptotic 
behavior of the diluted system. ._. 

In Ref. it was argued that the previous dataO for 
the cluster magnetization for system sizes L up to 48 are 
also consistent with the quantum criticality scenario. On 
a log-log plot, the last few points were fitted to a straight 
line, and the same exponent as that previously extracted 
from the staggered structure factor was obtained. The 
use of only a few system sizes in such a scaling is dan- 
gerous, however. It jjjeglects the slow curvature that is 
evident in the data.Ei With one more system size now 
available (L = 64), as well as increased precision for 
smaller L, the failure of the quantum critical scaling can 
be demonstrated even more clearly. Fig. ^l] shows a log- 
log plot of the cluster magnetization along with a line 
with slope-0.52, which was previously argued to describe 
the data.Ei The L = 64 point does not fall on the line, 
and also the smaller systems show deviations beyond the 
statistical errors. A slow upward curvature as L — > oo is 
evident. With no indication of a vanishing cluster mag- 
netization on the linear scale in Fig. [l(], the log-log plot 
is clearly not suitable for analyzing the data. 

In earlier -Monte Carlo studies of disordered Heisen- 
berg modelsali!3 it was concluded that the antiferromag- 
netic order vanishes in a quantum phase transition before 
the classical percolation threshold. With the current re- 
sults for much larger lattices at hand, it is clear that 
state-of-the-art simulations at earlier times were not able 
to reach sufficiently large system sizes for observing the 
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FIG. 21. The cluster magnetization vs inverse system size 
on a log-log scale. The statistical errors are smaller than th« 
circles. The line has slope 0.52, which was previously arguedEj 
to be the quantum critical scaling exponent. 



true asymptotic behavior. The essentially linear extrap- 
olations used to extract the critical point were therefore 
misleading. Similar work on the disordered half-filled 
Hubbard model on small lattices also have indicated that 
quantum fluctuations destroy the order before the perco- 
lation point. e3 In light of the behavior of the Heisenberg 
model, it would be interesting to repeat these calcula- 
tions using larger lattice sizes. For the Hubbard model 
it is currently not possible to reach system sizes as large 
as for the Heisenberg model, however. 

Experiments on quasi-2D cuprate antiferromagnets 
doped with non-magnetic impurities have in the past 
been able to reach only dilution fractions < 20%. til 
The doping dependence in .this region is in reasonable 
agreement with calculations a Recently, improved sample 
preparation techniques, involving simultaneous doping 
with Zn and Mg, have enabled studies of La 2 Cu04 all the 
way to the percolation threshold.Ej Neutron scattering 
measurements of the temperature dependence of the cor- 
relation length and the sublattice magnetization indicate 
that the order persists until p « p* f3 in accord with the 
behavior of the Heisenberg model discussed here. It can 
be noted that both the sublattice magnetisation and the 



spin stiffness extracted in the experimental agree reason- 
ably well with the curves extracted here (Figs. [l8| and |20|) 
at weak dilution but fall significantly faster to zero as the 
percolation threshold is approached. This is an indica- 
tion of additional quantum fluctuation mechanisms weak- 
ening (but not destroying) the long-range order, with 
likely candidates being frustrating nextjjiearest-neighbor 
interactions and 4-spin ring-exchange O The effects of 
these interactions are likely more pronounced in the di- 
luted systems. Random lattice distortions causing fluctu- 
ations in the nearest-neighbor couplings could also play 
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a role. 

In a system exhibiting a quantum phase transition as 
a function of some model parameter, e.g., the Heisenberg 
bilayerji^ certain types of dilution can drive an order- 
disorder transition before the classical percolation thresh- 
old. If the dilution leads to magnetic moment formation 
the phase transition is destroyed, however, since the mp=. 
ments interact and order even in the gapped phase.Ej 
In the bilayer, dilution of inter-layer dimers (two adja- 
cent spins on opposite layers) does not lead to moment 
formation, and a quantum phase transition can occur be- 
fore the percolation threshold (as a function of the dilu- 
tion fraction or the inter-layer coupling). A multi-critical 
point where the percolating cluster is quantuaucritical 
has recently been demonstrated in this system.Ej Quan- 
tum disordered ground states have also been found in 
Heisenberg antiferromagnets on nop-random fractal lat- 
tices, such as the Sierpihski gasket E3 
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